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1.1 Sources of Convergence Difficulties and Mitigation Methods 
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acceleration of the flow to sonic and supersonic velocities: typically, the Reynolds numbers 
are high and heat release terms are small. Associated with the various physical processes 
are specific time scales which describe the rate of evolution of these processes. In general, 
the time scales of the different physical processes are vastly different. This disparity presents 
a formidable problem to numerical algorithms as convergence of the flowfield solutions is 
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increases, this stiffness could be alleviated simply by increasing the CFL numbers. Figure 
1.2(b) shows the amplification factors for CFL u + c =500 (CFL U =5, CFL u _ c =490), where the 
stiffness no longer exists. This feature of improved damping with increasing CFL is often 
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Fig. 1.1. Magnitudes of the eigenvalues of the amplification matrix for the Euler implicit. 
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Chapter 2 



and the viscous matrices are given by. 
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acoustic eigenvalues revert to their standard values of u±c. In other words, when t = 1, the 
preconditioned system is essentially the same as the standardmon-preconditioned system, 
although some differences remain in the form of the time-derivatives. 
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standard Runge-Kutta algorithm only in the inclusion of the preconditioning matrix 
in the Jacobian matrix D. 



Stability Results of Preconditioned System 
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Fig. 2.1. Contours of the maximum amplification factor for the central-difference ADI 
scheme applied to the 2-D Euler equations. M=0.5. CFL=5, v/u= 1. 
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Chapter 3 
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contain a double superscript **. The second corrector (which is the final value at the new 
time step) is denoted by a superscript n + 1. The value at the previous time step is denoted 
by superscript n. 

In the PISO algorithm, the inviscid flux vectors are split into two parts, 
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with a similar expression for Bi. 

In keeping with the linearizations used in the PISO scheme, the matrices for the 
nonlinear flux vectors are (non-unique) multiplicative relations such that Em = AmQv and 
F n = BmQv The matrix Am is given by 


s in the linear flux vectors are centrally differenced, 
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3.2.1 Uncoupling the PISO Algorithm 
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This gives a diagonally dominant form of the equations of motion which is more stable 
than is the original form. 


3.3 PISO Stability Results 
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. 3.7. Contours of the maximum amplification factor for compressible PISO algorithm 
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rates at low particle Reynolds/low acoustic Reynolds numbers. The present preconditioning 
strategy, which involves scaling the time-derivative of pressure in the continuity equation 
shows promise in this respect because it provides a means of controlling either the acoustic 
speeds or the time-scale of the continuity equation whichever is relevant. We next consider 
the characteristics of the preconditioned system under viscous-dominated conditions. 
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We have discussed the characteristics of the Navier-Stokes equations in a previous 
section. At low particle and high acoustic Reynolds numbers, the dominant processes are 
viscous damping and acoustic wave propagation. The time scales associated with these 
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2y/2 while the limiting VNN is around 0.25. For the implicit ADI scheme, on the other 
hand, the optimal CFL number is between 3 and 5, while the optimal VNN is between 10 
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Reynolds no. limits o,*i /k k u v 3 /k condition # 



Table 4.2. Preconditioned Navier-Stokes, 
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Table 4.4. Preconditioned Navier-Stokes, e=Max(ej ; 
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first consider the Euler equations and treat the convective terms alone. We then turn to the 
complete Navier-Stokes equations where both convective and diffusive effects are present. 
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time step mav be obtained for all aspect ratios by specifying the minimum CFL 



time-step definition may then be enforced by rewriting the standard time-step definition in 
Eqn. (5.1) to select the maximum At in the two directions rather than the minimum: 

A( = M a x[^^, £™AL} (5.2) 
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therefore, are quite popular. MOCprocedures, though more complicated, are more rigorous 
since they draw directly from the mathematical character of the unsteady governing equa- 
tions. Specifically, the characteristics procedure dictates the precise number of conditions 
to be specified on each boundary, and provides the additional information needed at the 
boundary through integration of the appropriate characteristic equations. 



Past experience with applying various boundary condition procedures has indicated 
that all of the boundary condition methods (extrapolation and characteristic methods) work 
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Multi-stage schemes for the Navier-Stokes equations take the form, 



condition implementation and time step selection thus provides no improvements in Runge- 
Kutta systems. We do, however, note that the arguments are different for the Navier-Stokes 
equations and that proper viscous preconditioning and time step definition offers some pos- 
sibility of improvements at high aspect ratio. 
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number of the flow is 1 x 1(T 5 and the streamwise cell Reynolds number, Re& x , is 100. The 
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code described here as well as using an explicit Runge-Kutta algorithm are shown. The 
experimental values are obtained from the curve-fitted plot of the experimental data in Ref. 
[22]. As can be seen, the agreement between predictions and measurements is excellent. The 
slight discrepancy at Re = 400 arises because of the three-dimensional effects in the flowfield 
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. 5.6. Convergence to machine accuracy for inviscid duct flow with aspect ratio (AR) 
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Fig. 5.13. Predicted primary recirculation zone lengths compared with experimental values 
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The vector Q = (p, pu, pv,e, pYi) T , the standard conservative set of dependent variables. 
The viscous operator L is given as: 
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Fig. 6.9. Comparison of predicted velocity profiles with the experimentally determined 
values of Ref. [29] at various downstream locations for the multi-species, turbulent shear 
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Fig. 6.10. Comparison of predicted turbulent kinetic energy profiles with the experimen- 
tally determined values of Ref. [29] at various downstream locations for the multi-species, 
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Fig. 6.14. Comparison of predicted turbulent kinetic energy profiles with the experimentally 
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7.1 Formulation of Upwind Schemes 

7.1.1 Flux-Difference Splitting Procedure 

The convective fluxes are upwinded by standard flux- difference splitting methods. 
The numerical representation of the convective fluxes at the cell faces is defined by: 
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algebra, these schemes will be presented only for the Euler equations. Extension to the 

Navier-Stokes equations is straightforward. 

The two-sweep point Gauss-Seidel algorithm gives rise to a lower triangular matrix 
on the first sweep, and an upper triangular matrix on the second, thus giving an LU approx- 
imately factored system. 
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where S is as previously defined. 



Variations of the LGS algorithm include a single forward/backward sweep for each 
coordinate direction or multi-sweep schemes where several sweeps are performed at each 
time-level. Other variations of the four-sweep scheme may also be obtained by re-ordering 
the sweeps. A two-sweep LGS scheme is obtained by dropping Eqs. (7.5c) and ( 7 . 5 d). 

The upwind ADI algorithm takes the form of an x-sweep followed by a y-sweep, 
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terms in the LU system (Eq. 7.9) come close to cancelling. Note, however, that even for the 
direct inversion case we only approach a Newton iteration when the LHS and RHS employ 
consistent spatial differencing. 
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7.2.2 Convergence Results — Euler Equations 

To check the Euler stability results, a series of inviscid straight duct flows is computed 
with each of the upwind algorithms. The uniform grid consists of 41 x41 points and the grid 



aspect ratio is varied by reducing the width of the duct while keeping the length constant. 
This produces a uniform grid aspect ratio throughout the computational domain. For all 
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7.2.3 Navier-Stokes Stability Analysis 

For Navier-Stokes computations, convergence is influenced by both convective and 
diffusive effects. The overall importance of the diffusive terms for a particular computational 
grid is characterized by a non-dimensional viscous time step, the von Neumann number 
(VNN = i/Af / Ax 2 ). The relative importance of convection and diffusion is given by the ratio, 


d a 43 v 5 
J3 O o 

..1 3 Ii 3 

b> 5 S §3 v 

f|D Qh n ^ CO 

| §3. 1 f 

3 J: b g ? 

CD H3 «3 g ID 

> 2 8 2 S 
5 ! S « -S - 

v — ( w <D -+-» O 
a, o; <d c+h •- 

-S ^ * 0 g 

^ a g 'S ^ 

a |5 S-S 

ilaSs 

,T) A jS > o 

jj ^ g o 8 
Br & e ^ 

O £ -O <D 

o -C -2 -O 
_r "S S h 
# 8 b 3 s 
S d 8. fl a 


■ a 2 3 
=3 >, f -1 «j 2 

, «3 £ . Oh 3 

B (3 & — d 

co i — 1 O Q) J£* * 

<0 13 « u 2 ' 

co O 0) O ” 

5 ; _ <D bO d 

•“O _S F-O |^H r 

H § 2 - A 

8 g c * 

1 ° » s 1 ■ 

^ u2 3“ 

CD .2 O g *S t 

rt t-H o 

« <D & 

^ 3 pH £ 3 ' 

8 s g . 

e 3 to U J 1 

03 o *3 J2 d 

3 § & < 

a; ^ 

bO CO 5—* CO 

a? d -g 2 ^3 

eg I SS^- 

a? 

1 I-h <u ; 

s; S 

> +J CC 


CO 

o 




g g £ *s 

.2 8 8 g ^ 

s .a s> sz 

“ss-s 


03 03 

3 6 


03 

o 

d 

CD 

co 

CD 

fa 

Qh 

<D 

43 


I 

cd 

bJD 

g 

■a 

CD 

CD 

fa 

o 

fa 

cS 

d 

CD 

> 

w 



a <d 

d 

JJ s « . ^ . . 

5 ® u I .3 ^ O 

9 3 a 2 1 § 4 

£ g ^ g g 

g, ►» -§ g jj § a 

J <2 s. f? 2 * g 

4 * J5 * i* J I 

■ ■ ^ ^ CD 


* 8 
£ « 


(D (D 
CD N 


Cl & 

fa 


d 

0) 

eft , r~} 

3 So ^ *> +* 

■" , d \3 ^ d 

d CD «H 0 
O O cb 


O 03 


q 6 


<43 g 

' fa 


° 

CO CD *d 
- CD bo « 

O « -^3 

^ CO 

2 § ‘C '. 

.d ^ o (D fa o SJ 

Cd > ’^ g 

§ -§ a g * 8 2 

u ^ <+-t o 13 CJ 

o -2 U ^ OT w w 

‘ d o> 


8 a I 

3 ° 


d 

O 

3 


co ^ 


co 

fcg 

d 

CD 

fa 

d 

Oh 

CO 

d3 

fa 

43 

CD 

<b 

O 

fa 

Oh 

& 

fa 

CD 

> 

CD 

d 


CD 

in 

*> 

.s 

<D 

43 


sr s £ g 

.2 * w h 
“ _^J 2 
cd co rS 

W 4-t I X . d 

I - S - "1 
"■ § "E 1 §• § 

E W g 


CD > 
fa CD 


<D 


O 

CD 

CD 

< 


•fa d 


4j 


S CD 43 fa o ^ 

as 2 5 S ^ 2 § o> 

4h CC3 -H jH w 

S &, I 9 “ 
° -g a g 5 g 

— 1 Q «' -» 


t— M t-H 

CD rb 


<b ’ 
Oh 


bO co 

co CD 
cb 43 


fa 

.2 

03 

d 

O 

fa 

bO 

.8 

t-H 

* CO 

CO 

cb 


fa 

d 

fa 

d 

l 

.8 

fa 

03 

CO 

03 

d 

O 

*3 

.2 

449 

CO 

03 

d 

0 

& 

CO 

d 

0 

• 

♦ H 

JD 

13 

^03 

"cb 

03 

43 

d 


d 


0 

CD 

0 

42 

fa 

d 

d 

d 

U 


CO 

CO 


42 cr 
o w 

fa v , 

a 

S3 <3 
3 

8 .2 
CO -4-3 

*> d 

- § 

03 
fa 
Oh 

CO 

d 

o 

CD 
CO 


CO fa 

.S -2 


*n3 


J2 'g ^3 

13 2 « 

_ B ^ £ s 

.3 « 8 5P & 


03 


CO 


cb 

C-H 

4-3 

o 

03 

a 


4 s s a» !• 

5 03 rj .H O 

s § | -s « 

X g g, 3 A- 

g ^ ” 8 

E 8 .2 3 


a d a 

.2 .2 § 3 

o *53 

CD C/2 


I 3 


««>»<» 
- — 

bO 43 iz; 

h +3 rz; 


43 ■ 

bO 


f-i 

<2 

03 

6 

03 

43 

V 

co 


z 

z 

> 

x 

cc3 


£ 
3 

•5 M H a 

3 ° ^ • -a 

8 “ ^ S 

3 2 a 5 

^3 ^ a .3 

> ^ - 1 


o 

CD 

03 


*h 43 

c3 


S S 

8 % 

s o 

Cu I 

1.3 
I a 

Oh 


d 
o 

cD /-i »— i 

03 B P 

d ^ 
o- 5 p 


P □ .rH 

-38 
3 8^1 


■dI^ 

^ ” ^5 1 
<D 03 bC O 
d 43 *d 

■ ^ ^ 03 

f-H 

d 


03 .05 
43 +* 
*Z3 Cl 


<D 

£ 

CD 

43 

u 

CD 

13 

# S 

cn 

i 

w 

CO 

d 

d 

0 

d 

• *-H 

o 

Oh 


CO 

ci 


^03 U 

03 .3 

°? 6 
co 

« 03 

§ -5 

^ bC 

a -S 


ro 

i 


ft 

<1 

<13 


O 

Oh 

03 

43 


yj - O 

> >, -d 

^ g S 

2 o ° 
Dh -c 3 Q.O 
n ^ a # co 

O 03 V 

Irt Cb ^ 
(3 V o 
M ^ s vd 

t3 -2 OJ 

8 a 3 

rr ' (j 4-3 
« CO 

03 >, 

43 Zd 

S.i 

o 

o 


03 ‘-M d 

O o 03 

o .2 ^ d 

^ >> co .a 

Sl-s 9 


03 

43 

4-S 

03 

£ 

43 

o 

,p § 


03 cc3 H 

Ph fv* .43 


X 

• rH 

f-4 

d 

3 

bO 

d 

* 

d 

o 


'd 

d 

o 

CD 

03 


O 4 in 


cb 


42 

cc3 


B ^ 

ct O 


Oh 

03 

43 


to 

a 


co 03 »4H 
03 

CO 
0 
cc3 


^ 03 

U) 

o ^ 3 

— I X OT 


CO 


|-S^ 

JH 

M w « 
03 S j3 

£ J « 

rQ >1 
CO 


CO 


CO 


h g a> 

O^cS T3 
•<i2 2 


co M 

^ fa 

W Ch 
03 .d 
M 


4D 

cb 

4-3 

CD 


d 

£ 

o 

43 


d 
d 

CO 
03 

8 

•fi 8 

^ D 


03 ^ O 

D 

O d ^ 

II § *s 

CH ^ ,-s 

< 2 ^ 8 | 

^ v S® d 03 

03 x •"* Q. 43 

^ || ^ 03 m 

K II 03 r j^ * — ■) 

Oh ^ h4 

E§ ^ 
^ ^ q 4 fa 

lO w «3 4J 4-J 


o 

M L, 

Oh S 


CO 

fa< 

o3 

03 

Oh 

& 


d 

y3 

03 

"d 

03 

43 


co 

d 

o 

CD 

co 


3 £ 

' s 3 

tB 

co «3 id 

3 .2 ^ ^ 

v -g « 2 

r— < d 5> 43 

fat Q3 4-3 

O ^ * - 


"d 

03 

JS 

13 

fa< 


CO 

d 

43 


^ 03 ^ 

3 » ^ H 

u «3 H 

•r ^5 .2 ^3 
> hCl D 

-S fa g s 

o <+H cr* 72 

-d O 03 co 


•d 

d 

o3 

03 

fa 

d 

d 

03 

CD 

o 

fa 

Oh 


^ O | 03 

8 ”-§1 a 

'T H 0) y 

cn fa fa « C3- 

O PH V 4h 

o fa-H ^ CS 

>— I ^ o O o 

03 43 ^ d 

•5 -1 s ^‘S 

d o 


^*5 .3 3 
3 h 2 3 


I— < r—-H 


CO 

2 : S cc 

3 .> « 

H-j d 


o ,S co 


fat 

«3 

2 .-e &. -o ^ 

® "O o c o 

■2 § -d g 1 

^ y U % 1 


CD 03 

H fa a ^ fa 

<3 p 


CO 

d 


CO 

3 


o 

*3 


d z 

CO ^ 

03 . 

43 pj 


<1 ^ 

S 

4 O r 1 

o d • y 4^ S 

fat c; 33 cj ^ 

Oh S3 d s 

2 O g d fe 

do 
o S 

y n 


03 

l 

_ 43 

d u 
d co 

« r- d 

sgjl B ° 


•g ^ CO - 

d - d -5 
03 o o d 
03 ^ y S 

^ ii ■ a 

H > 


> 

o 

bC^. 

d ^ 
• d rt 

Oh 3 

S 


_ > 
O 

03 d 
d °3 


co 

*T 03 
.2 43 


fa 

o 


8 5 9 


W *• -77 

Is 

.2 ^ a *g 

rt «3 d 

M s SP 03 

4-3 0 


CO 

^ o 

>, 


03 

£ 

03 

43 

CJ 

cn 

'aJ 

d 

’S 

cn 

i 

CO 

xn 

d 

cd 

0 

03 

d 


CS 

CO 

CS 


T 5 Oh d co 

O ^ Q 3 »-H 4J 

£ 03 fa .3 

sr « 


Oh B 

W d 
^ to 

43 33 
bO ^ 
’d co 


03 


co ^2 03 

fa d fa 

° o ^ 

■p c <5 

CD d 4f 

cb ^ d 

fa - ^ .S 3 

o 03 s X ^ 

■w CO O (3 +■> 

qj •»-» w co 

w fa ^ 

4 J 4 h aj jS 1 03 

"d H y 03 43 

8 . & d ^ 

2 3 a ^"S 

| 3 ij 

<H ^ fa ^ 

S .SP 2 tK "3 


Z5 2 «3 fa4 
CO > d 
d II ^ 43 

! L 03 o 

d h ^3 w 

O t-4 4J 03 d 

O tfa 43 03 

I’S ^ r il 

« , fa ^ > .2 

”T 0) 0 a> ex 

H ^.2 6^ 

3 u o -g « „ 

S.2 2 «- -6 •" 

fa 43 03 M ,2 

g 1 2 cS w 

■+ 3 c2 +» c . 

•3 2 0 ^ lz 

3 » 2 z 

^ d 3 ^ 

y d ^ rT 

5 03 fa 

Oh Q_, q fa - h 

. w 0} tS 

03 ^ 0; ^ 03 
g d ^ d d 
- 10 *3 

4-3 


d 

fa 


d . 

Oh 

” bb 


45 ... 

«3 [i P fa 


03 

43 

H 


a 2 
| II « 

w ^ 

2 z -3 

3 " 


CO 03 

>> & 

fa 


^ d b 

2^4 . r% 

H 2 CO 0 

c -S g w .2 


3 

d 

8 


Oh 

o 


03 

> 


^ w d d (u 5i 

&> ^ ^ *: 8 *5 

s a 


d 

93 

fa g 03 

x 4 d 


bO 

S 

44 

CD 

CO 

Id 

H 


4 4^ 03 

<u g - 3 

g 3 1 a 1 ^ 

fa d H fa 
-fa d S O 
cj .d d 13 

0) 4-5 fa fa 

P, Oh d 03 
co O 0 0 


cb 

fa 


d 

.w ^ O d 

Oh Oh CD h4 


o 



= g'|>l g 

a; g od *3 .o m 

d « ^ *£ ^ cn a 

” t£i y -a ^ ^ ^ g 
►>. - d d d ^ .d 

_j d bO u |ij ^ h j 
dMCDEoddpp+d 

•Jig 'in 8 

■« 8 3 ' ° ^ g 

§ £ * 8 ^ g> 

t 1 • 3 S m H ->j m jrj 

3 r3 -E -r <t _ a; <u 


■i - 3 r - 1 X 

3 M ° 


t3 

o 

CJ 

co 

d 

. ^-1 

>3 

pp 

pp 

d 

o 

o 

CO 

<u 

u 

d 

’’d 

d 


d 

J2 

pp 

CO 

d 

> 

c/T 

pp 

CO 

.S 

co 

d 

d 


<u 

0> 

^3 

fai 

n 


*■■ co 

<D XJ *a O H 3 

d D r^ o 0 ) 

g § H N “ 

od ^ 

o V (M cd 

CO b Oh »h 

Co . m d 


P Tj D 

m ^ S 

d d .d 


O ^ CS 1 CD _ ^ 

CO H I Oh (h O’— < ►_,, 

W » N m« tO *ZL 

a ? 7 o 1 > 


„ « s a s?; x ai 

5 S 

»-i' g, la 3 g j 
B H • SP «3 « fa 

m CU “ O 


g H . £° <3 

JK So S) J2 . “ 

'o (3 £ (h H <l> 

S’ - " _ <3 fa _q 
ft 3 ,g O - 

. <D <u ^ 

^ ri w > ^ — 

<D ^ +* O £ 

d «- d n . H 


3 .S - 
S E £ 


cp _Cj trj r> p— » _j p/ ^ 

<D *> q g 5 ^ ^ 

P! >73 ^ b 0 O pp pp 

I f S a | j ’ 

£&<»•- dfag? 

-£ 

7 O ^ pp d - *- 3 co 

d P «3 PH ^ < rO 

o A ^ o ^ p d 

« ® »i« 3 - • E 

•gera ^ e 8 g' 

8 o -a » >. x ^ 2 : 

ED .3 2 Sa ^ © T3 

J I 7 8 ^ 03 . 

(o B — o -js o 

_c S •> 2 o <3 p 

pp (_< J> G p (p co ■ 


$ >,3 -5 £ 3 

&>*?<§•§ 
4 J £ (n SS ^ 6 

CO 0d ^ ^ g 

^ pp CD 9 k n 

•^ 0 S5o « 

^ « p % > £ 

A 0 ) > .2 D s 

D ^3 a"? - 

P 7 P 8 >, 8 

^ 8 8 ^ 3 E§ 

o .s B « *2 S 

^ W gg ft) o 

B . i B J c 

E 2 x Z o 

8 2 7 -5 S 2 

8 *6 E ^ .2 ^ 

fl3 CO (p d) P3 0 

rt l-H & ^ g <U 

g q ^3 ^ a 

^ J-j a m-, g g 

^ «3 ^ a3 ^ 

W5 g U +5 a *k 

o ^ | 3 .o & 

r-r-y ] QJ < ^ > i-^H 

1 $ ^ ^ 2 ^ 7 

0 ^ ^ a < 

m O a bp.S 

.S j ri d o * 

^ ,.S a° 

^ a) ^ “ vh 

o 0) «3 <L) II 

*§ - .a a g h 

tp p u ° <f 1 

t3 ° «S w 

»p Qj O O fyJ 

"§ 2 d >,-a? 


.a o 8 

W 3 P ^ 

° g j> 
l-> P 3 > 

» *» ^ 
O M 
• p o • 

pp ft 4 P CO 

q; 

^ CO m 


» g -s i’ a -g 

I g’.rsi! 

.52 > P ‘h "O o 
> S W <3 C 5 o 

^ o g ^ u o 

J ^ p d — 2 C0 

Oh P CO • tP d •- 

wig 8 
g co *c ^ a ^ 
£ a> &> p * s £3 
^ p o -2 § 

i?^ ^SB « 
3 *« v, ^ • "6 


(jp >> 


S 4 a 9 ^5 
p ^ -s i§ s 

3 « g “ e 

cc 3 Oh O 

w 73 S 8 H 


d a) 
^ CO cn 

IS S 5 

o 

. |3 ^ 

i^li 

^ <L> o O 

II E 12 - 

Cd O m o 
p «3 P 
< o ^ ^ 

CO ^ pp 

P Z m 

T'ob ^ m 
o O > g 
d p p * 

^■5 o z 

b *+- E> 
II , o o 
It ^ _ <u 

pS co .2 -O 
7 ? ^ p ^ 

^ 3 2 -a 

§ 

> s 8 Js 

a &U3 

O fp o 

,»P 0) P 

^ > CD % 

co d r-C co 
(D n p 

CO 3 .. CO 

d w O «3 

OJ co O > 

H ^ PS ^ 

O > *<p _ 

dj ocnJ 

d ^ (jp 

co . 
o O 
< 1^2 i 

c^ o- 7 e 

m o M 

3 3 7 2 

P bp ^ "5 

0 E -| ^ 

d ^ d 


(PDiOTrjcOPP'd^ 

is a g | g g | s -5 

1*3 s-s| >. 

ilESIJa *g^ | 

° 5.5 s i H “ 

S ^ d ^ ^>7 >" u.S‘CO 

S c fl «> e 2 r 

s s-sil 


-2 ell a « *| is 

S3 § t 3 d i a -2 ^ ^ l : ^ 1 | 

>pg ^^'-. d ’5 . mC S 5 

^S'Sp^ a d 

-o'! > 8 « .3 8^ B 


s < -> CU cp. ^ cy pj co 

fcsi|lS|S-£8g| 

u 5 'i 1 Is 3 a* §■ 

O- 3 o 3 B. bC 5 no 
^ 7 ^ 0 G C m fe— 2^0 
^ fa o •Op-->Sa>f 1 

o n V * u ^ sS'goHB 

<«^.s s 8 a 

OOcjCP M ^ U OT) V p 

"O ^ O 3 fc 

$ •■§ 5 p ts | :s s § ^ 8 ) g sa 

^ IS J 8^5^ I 

^d S rd -2 * Pp ^ xo QJ co rP^ U 

d ^ ^d co 52 ^ S 4 g *>> o c^ ^ ^ 

2 Sp 53 . 52 °g“|po£ 

J 8 I'^““p cc3 § -s og-S 


fa( t>^ ^ co 

2 g3 I 

25 0 v Jtt 

Z U J 3’5 

> £ ~ ” 

■ — ^ 0) 

j 2 S ^ 

fa P O O 

fa> -g * s 

[ ph d 7 p 

° |f> E ^ •§ 

.2 0^7 ^6 

pp 0J t— 

•-"$£«! 

1 § ^ a.g 

g p o •- o 

S £ ~ 43 

131.9 2 

o m .a a p 
nj P 0 S *2 
_S 7f cr g o 
■5 > oj o « 
-- ^ *w 0 > 

92 « m g 

E 2 g .2 S 


E 2 s 
« S .2 
•f? p ^ 


” 2 ^ 

CO d ^ 

d co « <u 

2 ■" p 3 

d Si co co 


p p d co d JrU p fafa z v, 

2^2 co.9^W)5^^ *h On 

^2^a>Ii ^js-2s« 

2 S^3 3 » ° S “ I? « i 

e sp- 2 r: , s>-- « cS S 5 : 1 :i 


p p ^ 

p § ^ 
E B 'Z 


p 

d 

^d 

pp 

pp 

d 

CD 

O 

<D 

Jd 

p 

> 

r-d 

fp 

o 

o 

pp 

o 

Up 

co 


oo 

CO 

pp 

CD 

o 

CO 

CD 

6 

pp 

pp 

<D 

d 

<D 

f N 

pp 

O 

W 

</> 

Oh 

d 

CD 

CD 

f— H 


bCQ 

3 < 

<D 

02 


o 2 «3 s 
m O g 8 

p g p ° 

pj d , 22 


2 i ^ 

■g i p -§ 

g E *= 3 

|5S J1 


2 , §’ s 's.S| 

O ^ Oh QJ 

S *7 J 

§ s S p Q & 

« J I 

5 1*^1 8 

5 -£ < *= 

77 ^ co 

^ co d 1 ^ 

1-0 rn CO CT 3 d 


convergence as the aspect ratio is increased. 



7.2.5 Timing Study 



CODE | SCALAR | VECTOR 
ADI I I ^ 1 

LGS-2 1.08 1.15 

LGS-4 1.63 1.54 

LU 1.08 1.16 


Table 7.1. Relative computation times for the ADI, LGS-2, LGS-4, and LU schemes^. 
Times referenced to the ADI scheme: 4.0 x 10" 4 sec/step/grid point (IBM 370), 4.8 x 10 5 
sec/step/grid point (CRAY Y-MP) 
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amplification matrix for the Euler equa- 
tions: Direct Inversion, (I/III), M=0.7, 

flow angle=45°, AR=100, CFL X =1000, 
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Fig. 7.5. Maximum eigenvalue of the ampli- 
fication matrix for the Euler equations: LU, 
(I/III), M=0.7, flow angle=45°, AR=100, 
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Fig. 7.8. Convergence results for invis- 
cid straight duct flow at various AR’s us- 
ing the LGS-4 scheme with I/III differenc- 
ing. M=0.001, min-CFL=20. 
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amplification matrix for the Navier-Stokes 
equations with min-CFL/max-YNN pre- 
conditioner: LU, (I/III), M=0.001, flow 

angle=45°, AR=100. Re& x =\0, CFL :r =10 6 , 
VNN„=10 5 
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Fig. 7.12. Maximum eigenvalue of the 
amplification matrix for the Navier-Stokes 
equations with min-CFL/min-VNN precon- 
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PRECONDITIONING FOR UNSTEADY COMPUTATIONS 
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We use von Neumann stability theory to enhance our understanding of these issues and^to 
guide the selection of the preconditioning matrix. We then use various test cases to confirm 
stability findings and to assess the capabilities of the dual-time scheme. We note that while 
our focus is on implicit algorithms, the formulation of the dual-time scheme is applicable to 
explicit algorithms as well. 
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significantly within the solution domain because of local grid-stretching. In the following 
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A more general procedure for obtaining iterative solutions to unsteady flowfields, and 
one which contains the approximate Newton method as a subset, is the dual-time approach, 



wherein a ‘pseudo’ time-derivative is introduced in addition to the physical time-derivative 
in Eqn. (8.1): 



where again central-differenced ADI is shown for simplicity. Note that Eqn. (8.4) has been 
derived by linearizing the governing equations is terms of the new dependent vector Q v . 
Also, S is given by: 
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8.2.4 Analysis of the Dual-Time Scheme 

Von Neumann stability analysis of the one-dimensional version of Eqn. (8.3) is useful 
in enhancing understanding of the dual-time scheme because the amplification factors of the 
scheme may be derived analytically. The amplification matrix may be written as. 
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(steady preconditioning), one of the eigenvalues approaches unity as Mach number becomes 
low. This will correspond to poor convergence rates in the inner iterations. These conclusions 
are in agreement with the stability findings in Figs. 8.2 and 8.3. 



Alternately, in the limit of large physical CFL s (the limit of steady-state), the eigen 
values become: 
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LGS, this results in poor convergence rates. However, when e = M (steady preconditionmj 
the three eigenvalues are of the same order of magnitude and the difficulty is resolved. 

Examination of the eigenvalues in Eqn. (8.6) suggests the following definition of e 
as to maintain them of the same order of magnitude under all conditions: 


The above expression depends on the wavenumber through the term S x - Since in most 
problems, it is the low wavenumbers that are typically the hardest to damp out, we tune 
the e definition so that it is optimum for the smallest wavenumber in the domain, i.e., 
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in an optimum fashion— one grid line per physical time-step. At this choice of time-step, 
however, the approximate-Newton method becomes completely ineffective. In contrast, the 



8.3.3 Spatial Instability of Mixing Layer 

As a third and final demonstration, we consider spatially growing instabilities in a 
simple mixing layer. Figure 8.13 shows a schematic of the physical set-up. We consider a 
hyperbolic tangent shear layer whose mean profile is given by: 
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From Fig. 8.14, it is evident that the unsteady preconditioner (in Eqn. 8.7) performs 
?t. reaching 1 machine zero in about 600 iterations. The no-preconditioning case may 
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show good qualitative and quantitative agreement except near the exit boundary, where 
boundary reflections appear to impact the solution contours. Quantitative assessments of 
the accuracy of the predicted growth rates and frequencies are reported in Ref. [47]. 
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Fig. 8.2. Stability map for the approximate- Newton scheme (CD/ADI) for M — 0.001 
Lower Triangle: CFL u+c = 1. Upper Triangle: CFL U =1. 
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Fig. 8.3. Stability map for the dual-time scheme with steady preconditioning (CD/ADI) 
for M = 0.001 and CFL r = 3. Lower Triangle: CFL u + c = 1. Upper Triangle: CFL U = 1. 
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Fig. 8.6. Inner-iteration convergence shown for several consecutive physical time-levels 
for tracking acoustic oscillations in a cylindrical chamber. M = 0.3 and CFL u +c = 1- 
Convergence is shown both for the approximate-Newton scheme and the preconditioned 

dual-time scheme. 





Fig. 8.8. Initial specification of vorticitv for the propagation of Lamb vortex in a straight 
channel. 





Fig. 8.9. Convergence rates of inner iterations for a single physical time-level in the propaga- 
tion of Lamb vortex through a straight channel. M — 0.001. CFL u + c = 10 and CFL U = 0.01. 
Scheme is central-differenced ADI. 
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Fig. 8.11. Convergence rates of inner iterations for a single physical time-level in the 
propagation of Lamb vortex through a straight channel. M = 0.001, CFL u+c = 1000 and 
CFL U = 1.0. Scheme is central-differenced ADI. 
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Fig. 8.14. Convergence rates of the inner iterations for a single physical time-level for the 
computation of a low-speed unsteady mixing layer (CD /ADI). 










Fig. 8.16. Instantaneous U-velocity and V-velocity contours for unsteady mixing layer 0.01 
s after start of computation. 
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the smallest grid dimension. The CFL number associated with this direction is thus speciii^d 
at an optimum value, while the CFL numbers corresponding to the longer cell dimensions 
are forced to be very small. As a result, flow-field disturbances propagate slowly and are 
damped ineffectively in the longer cell direction, which results in poor convergence. In 2-D, 
the two factor ADI scheme, allows the local time step size to be specified on the basis of 
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The preconditioning matrix, I\ is the same as that which was previously given in this report. 

Applying Euler implicit discretization in time yields the following unfactored algo- 
rithm in delta form: 
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diagram looks well conditioned over most of the wavenumber space. However, it approaches 
unity for purely streamwise modes (a consequence of the extremely low CFLj), which is the 



classical problem with high aspect ratio grids. Thus, the three-factor ADI scheme gives very 
poor convergence for grid cells of Type II no matter how the time-step is chosen. 
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second pseudo-time derivative: 
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choice of Q° is taken. Thus, the stiffness associated with the streamwise modes, as weir&s 
the high y and z wavenumber modes, has been alleviated through simply making the correct 
choice of the starting vector for the dual time iterations. 
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aimed at understanding these two more complex systems is needed before uniform, reliable 
convergence can be expected in these types of grids. 
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Abstract 



incompressible Navier-Stokes equations. Test cases involve laminar and turbulent flows with 
severe gnd clustering and flow separation. The results are compared with those of a basic explicit 
Runge-Kutta (RK) time-stepping iterative algorithm and with the Implicit Residual Smoothing 
(IRS) and the Distributed Minimal Residual (DMR) acceleration techniques. The new acceleration 


scheme is shown to be superior to these methods especially on highly-clustered grids. 
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Here, we plan to elaborate on an optimum extrapolation technique in order to find optimum 
relaxation parameters for iterative algorithms. An extrapolation procedure based on the power 
method and the Minimal Residual Method (MRM) was successfully applied [4] to the multi-grid 
algorithm. The MRM uses equal optimal weights for the corrections to every equation in a system 
of partial differential equations that is solved, but it has not been shown to accelerate the scheme 
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where x is the time, Q is the solution vector, E and F are the inviscid flux vectors, H is the vector of 



Here, U and V are the contravariant velocity components, u and v are the Cartesian velocity 
components, p is the pressure, l* and r\ are the curvilinear, non-orthogonal, boundary conforming 

grid coordinates, J is the determinant of the Jacobian transformation matrix, J = det(d(|;n)/d(x,y)), 
and p is the user-specified artificial compressibility parameter [17]. The term d( jj^/dx was added 


to the mass conservation to make the entire system non-singular. Adding the non-physical term. 

does not affect the steady-state solution since all the time derivative terms vanish as the 
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Here, R is the original unsmoothed residual, R is the smoothed residual, 6 is the second order 



Distributed Minimal Residual (DMR) Method [5-8] 






b . b 





optimum co's to minimize the L-2 norm of the future residual of the system integrated over the entire 



domain. The present formulation of the DMR uses the same values of the N x M optimized 


o 







optimum amount of corrections necessary to minimize the overall future residual. The basic idea is 

that the residual at a grid point depends on the values of the solution vector Q at the neighboring 
grid points used in the local finite difference approximation. The rate at which the residual, R m (m 

= i M: number of equations in the system to be iteratively solved), changes with Qg (s = 
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Using the first two terms of a Taylor series expansion in the artificial time direction, the residual for 



Subsequently, the future residual at (t+n)+l can be approximated by 





162 








In equation ( 14), R's and a's are known from the past iteration results. Since each a is assumed to 






In the general case of a system (1) composed of M partial differential equations, there will be M 
constant optimum values of a. This means that the system (15) will become a full MxM matrix for 


M unknown optimum a's. Also, each summation in the system ( 15) which in this example is from 

m =1 to 3, will become a summation from m = 1 to M. 

Each application of the SBMR creates significant disturbances in the eigenvalues of the iterative 
matrix. Thus, the optimum values of a’s are determined and applied only periodically followed by a 
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constant over the whole computational domain, it should be classified as global SBMR method [12, 






o 






simultaneous system of equations (17) yields a block penta-diagonal matrix equation for (jmax- 
jmin+1) x 3 optimum a's where each block is a 3 x 3 matnx. In the general case of a two- 
dimensional system (1) having M partial differential equations, the block penta-diagonal system 
( 17) will have blocks of M x M size. In the case of a three-dimensional system (1), the system (17) 


becomes a block septa-diagonal matrix. 
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outperformed other acceleration schemes. When the maximum cell aspect ratio was increased to 


if, It" 




gradient along the inner part of the turning section resulted in an open-end now separation 
downstream of the turning section (Fig. 8b) which was successfully predicted because we used the 

correct form of the nonreflecting exit boundary condition. The LSBMR method used the results at 
30 iteration levels apart. The convergence histones for the maximum cell aspect ratio of AR^ - 



developed and applied to accelerate the convergence of the explicit RK algonthm for 



ible Navier-Stokes equations. The methods are easy to comprehend and to implement 
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Acceleration of Iterative Algorithms on Highly Clustered Grids 
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ACCELERATION OF ITERATIVE ALGORITHMS 



acceleration schemes 





and three-dimensional laminar incompressible flow Navier-Stokes equations, two-dimensional 
turbulent incompressible flow Navier-Stokes equations and two-dimensional compressible flow 
Euler equations. The new sensitivity-based acceleration methods demonstrated superior 
performance in all test cases that involved severe grid clustering and grid non-orthogonality and 


included laminar and turbulent flows with closed and open flow separation. 



effectively reducing low and high frequency errors yielding impressive convergence rates, its 



efficiency is significantly reduced on a highly clustered grid^ 7 and it is difficult to implement 
reliably 8 . The preconditioning methods 8 ' 10 , although very powerlul in alleviating the slow 
convergence associated with a stiff system for solving the low Mach number compressible flow 
equations, have not been shown to perform universally well on highly clustered grids. 

Optimum extrapolation techniques represent an entirely different type of convergence acceleration 
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order accurate central differencing throughout the study. Characteristic boundary conditions will be 
used at the inlet and non-reflecting boundary conditions will be used at the exit boundary of the 
flow domain. We will be seeking accurate steady-state solutions only. 



II. Global SBMR Method 



R with respect to each component of Q. This information is utilized to effectively extrapolate Q so 
as to minimize the future residual. The rate at which the residual vector components, R (r = 
l,...,rmax: number of equations in the system), change with the solution vector components, Q^, ,(m 

= 1 ... M: number of unknowns; s = 1,...,S: number of surrounding grid points directly involved in 
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assumed constant at every grid point for the global SBMR method, as 
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In equation (7), R's and a's are known from the preceding iteration levels. Since each a is assumed 

to have the same value over the entire computational domain, equation (7) gives a system of M 
simultaneous algebraic equations for M optimum ctp ct 2> • • * , ct^. 




For the general case of a system composed of M partial differential equations with M unknowns. 
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Navier-Stokes equations as an example. The system has three equations (rmax - 3) and three 

unknowns (M = 3), that are pressure, p, and the local velocity vector Cartesian components, u and v. 
The acceleration coefficients for those unknowns are a a u and a y respectively. If the grid lines 
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blocks [A], [B], [C], [D] and [E] are coefficients of a J ' 2 , a J ~ 1 , a-*, a J+ 1 and a J+2 , respectively. Each 
column of the blocks [A],[B],[C],[D] and [E] corresponds to a p , a u and a y , respectively. 

Application of the LSBMR method to a three-dimensional problem would result in a large. 
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compressibility method 26 . The Courant-Friedrichs-Levy number and von Neumann number had 
maximum possible values (CFL = 2.8 and a = 0.4) in all test cases discussed in this paper in order 

to achieve maximum performance of the basic RK scheme. 
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8). It is interesting to notice that TSS did not this time have any effect when applied to the basic 
Runge-Kutta scheme, while achieving impressive acceleration (Fig. 8) when applied together with 
the LSBMR method. To test the acceleration methods for situations with flow separation that 

stretches through the exit plane 22 23 we computed laminar flow in a U-shaped channel with 
computational grid symmetrically clustered towards the channel walls. With AR max = 200, 



LSBMR method offered the fastest convergence (Fig. 9) in this case, while IRS required more time 
than the RK scheme although coefficients of the IRS were chosen for maximum performance . 

Performance of the SBMR method was next evaluated against the DMR method and the RK 
scheme for inviscid compressible flows. Convergence rates of existing algorithms are very slow 
when integrating compressible flow equations at low Mach numbers^ because of the excessive 
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the basic RK iteration scheme became much faster than for the low subsonic case (Fig. 1 1) and no 
convergence acceleration was obtained by either the SBMR method or the DMR method. Since the 
present formulation of the SBMR method does not account for the artificial dissipation terms, the 
convergence of the basic RK scheme with smaller amount of artificial dissipation for the transonic 
flow becomes slower 2930 , while the SBMR method still offers convergence acceleration (Fig. 12). 



B. Three-Dimensional Flow 

To compare performances of the SBMR, PSBMR and APSBMR methods, a fully developed 
laminar incompressible flow was computed through a straight duct with unit square cross section 
having length of ten units. The grid was clustered symmetrically toward the solid walls and a 
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performed with the basic iterative algorithm that are involved in the evaluation of the change of the 
solution vector. In the case of two-dimensional incompressible viscous flows without severe 
pressure gradient, the SBMR and LSBMR methods significantly accelerate the convergence of 
iterative procedure on clustered grids with the LSBMR method becoming more efficient as grids 
are becoming highly clustered For a two-dimensional incompressible laminar flow, the SBMR and 



the LSBMR methods offer fast convergence regardless of the value of the artificial compressibility. 



successfully reduces the computational effort for solving a three-dimensional, laminar flow through 
a straight duct without flow separation. The current formulation of the sensitivity-based methods 
do not include the effects of the artificial dissipation due to the complexity of such formulations. 
The general formulation of the new acceleration methods is applicable to any iteration scheme 
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